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relatively  large  particles,  it  can  be  driven  across  the  Atlantic  Ocean  until  the  south,  central  and  north  part 
of  America. 

-  Optical  Parameters 

The  desert  dust  is  composed  by  highly  irregular  shaped  particles  that  cannot  be  described  as 
spheres.  Our  technique  is  able  to  describe  these  particles  as  clusters  of  spheres,  possibly  characterized 
by  different  refractive  indices,  whose  arrangement  is  chosen  so  as  to  simulate,  as  best  as  possible,  the 
real  components  of  the  aerosols.  This  approach  offers  many  advantages:  it  allows  us  both  to  determine 
the  correct  scattering  pattern  that,  even  after  the  random  averaging  procedure  over  the  orientations  is 
performed,  still  resembles  the  asymmetry  of  the  constituent  objects(1 ,2);  both  to  follow  the  change  in  the 
optical  properties  when,  due  to  natural  processes,  the  single  scatterer  changes  its  morphology  (3,4). 
Furthermore,  the  T-matrix  approach  enables  us  to  determine  analytically,  when  necessarily,  averages 
over  non-random  orientational  distributions,  once  an  appropriate  weight  function  is  introduced  (5,6). 

To  simulate  the  configurationally  anisotropy,  we  choose  different  cluster  models  composed  of  4  or 
5  spheres  arranged  as  in  Fig.  a. 


At  1=351  nm  we  used  for  the  refractive  index  n=1.5225+i0.006. 

-  Size  distributions 

The  desert  dust  ranges  over  several  orders  of  magnitude  in  size  from  10  pm  to  10  pm.  The 
observations  suggest  a  size  distribution  represented  by  means  of  a  trimodal  lognormal  function;  however 
the  third  mode  that  has  a  modal  radius  of  5  pm  can  be  neglected  as  it  represents  a  size  distribution  of 
particles  that  have  a  great  settling  rate.  For  this  reason  we  consider  a  bimodal  lognormal  distribution  with 
the  following  expression 
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In  this  formulation  Nj  is  the  total  number  of  particles  per  unit  volume,  n  and  g\  the  distribution  of  the 
geometric  radius  and  width,  respectively.  The  parameters  we  use  for  the  distribution  are  derived  from 
literature  (Table  2).  By  varing  the  distribution  parameters  in  a  random  way,  we  compute  5000  bimodal 
distribution  that  we  use  to  perform  the  integral  necessary  to  compute  all  the  possibly  a  and  b  values  for 
our  model.  In  fig  (b)  we  show  some  of  the  distributions  calculated  with  three  sets  of  the  parameters  (n,  cn  , 
Ni)  and  (r2,  a2  ,  N2)  chosen  in  a  random  way  in  the  range  indicated  in  table  2. 

Table  2.  Parameters  for  the  size  distribution  of  the  desert  dust  as  reported  by  different  authors. 
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-  Results  and  Discussion 

The  simulation  that  is  still  in  progress,  considers  aerosol  scatterers  modelled  as  cluster  of  identical 
spheres  whose  radia  are  chosen  so  that  the  radius  of  the  sphere  surrounding  the  whole  cluster  ranges 
from  0.02  pm  to  2  pm,  even  if,  with  the  acquisition  of  more  powerful  computers,  we  plan  to  extend  the 
range  up  to  5  pm 

We  present  here  the  preliminary  results  obtained  for  the  first  cluster  shown  in  fig  (a). 


Fig  (b)  Fig  (c) 

In  fig.  (c)  we  report  the  scatterplot  in  a  loglog  scale  of  the  extinction  coefficient  vs  the  backscattering 
coefficient  obtained  for  the  desert  aerosol.  It  can  be  observed  that  we  obtain  a  dispersion  of  the  values 
that  we  are  going  to  compare  with  the  simulation  that  consider  spherical  scatterers.  At  a  first  exam  these 
values  present  a  more  large  dispersion;  feature  that  seems  to  be  an  interesting  peculiarity  of  the  model 
with  which  it  would  be  possible  to  analyze  more  correctly  the  scattered  light  generated  by  the  interaction 
of  the  Lidar  radiation  with  the  atmospheric  constituents. 
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Neural  network  scheme  for  the  retrieval  of  size  and  index  of  refraction  of  microparticles  from  the 
angular  scattering  pattern. 

The  retrieval  of  the  structure  of  the  atmospheric  constituents  from  experimental  data  is  one  of  the 
most  interesting  problems  that  have  been  faced  in  these  last  years.  As  the  most  recent  instruments  are 
able  to  sample  the  atmosphere  giving  a  large  amount  of  information  whose  inversion  require  high 
performance  computers,  has  grown  in  recent  years  the  necessity  to  develop  new  retrieval  algorithms  that 
would  not  rely  on  detailed  knowledge  of  physical  processes  involved  and  fast  enough  to  be  implemented 
in  computers  with  low-cost  hardware.  Because  neural  networks  generally  exhibit  these  properties,  we 
implemented  this  algorithm  to  retrieve  the  size  and  the  index  of  refraction  of  a  dispersion  of  spherical 
microparticles  from  the  one  angular  scattering  pattern. 

In  the  framework  of  the  neural  network  a  scattering  problem  can  be  considered  as  follows:  the 
intensity  of  the  scattered  light  measured  at  an  angle  9  coming  from  a  single  particle  embedded  in  an 
external  medium,  is  a  function  of  r  and  m,  the  radius  and  the  refraction  index  of  the  particle  respectively 

l(9)=F(9,  r,n) 

Given  a  vector  x  containing  n  measured  values  of  I  for  n  different  9  angles  (these  represent  the  input 
data),  the  inverse  problem  consists  in  determining  an  inverse  function  f=F1  that  gives  us  the  information 
on  r  and  n  (the  output  data). 

-  The  training  data 

We  have  calculated  a  set  of  vectors  x  containing  the  51  elements  representing  the  intensities  obtained 
at  X=632.8  nm,  by  varying  the  parameters  r  and  m:  these  vectors  have  been  used  as  input  quantities  for 
training  the  network.  The  radius  of  the  particles  has  been  varied  in  the  range  of  50-500  nm  and  in  Table  I 
we  report  the  values  of  the  real  and  imaginary  part  of  m2  used  to  obtain  1(9). 

As  one  can  observe  in  Fig.1  the  input  data,  even  in  the  case  of  sphere  with  the  same  dielectric  constant, 
presents  a  large  spread  of  the  numeric  values  representing  the  intensities  (y  axis)  so  it  has  been 
necessary  to  perform  normalization  and  scaling  in  the  range  [0,1]  of  all  data;  moreover,  for  a  better 
performance  of  the  network,  the  inputs  calculated  in  the  angular  range  20°-170°  were  disordered.  In  Fig. 
la  we  plot  the  same  data  represented  in  figl  after  the  procedure  of  normalization  and  mixing. 

•  The  training  algorithm 

We  have  employed  a  feed-forward  neural  network,  also  called  perceptron.  The  type  of  training  algorithm 
employed  that  yield  the  best  results  is  based  on  the  Levenberg-Marquardt  method.  This  algorithm  gives 
fast  convergent  results  even  if  it  requires  a  greater  amount  of  computer  memory.  To  take  advantage  of 
the  neural  network’s  nonlinear  approximation  capabilities  requires,  in  addition  to  an  input  and  an  output 
layer,  one  hidden  or  intermediate  layer.  The  activation  functions  used  for  the  hidden  layer  are  the  tansig 
function  and  for  the  output  layer  the  purelin  function.  Moreover,  the  exact  size  of  the  hidden  layer  has  to 
be  determined  empirically. 

-  Conclusion 

As  shown  in  Fig  2,  the  most  successful  network  configuration  found  in  this  study  uses  51  inputs  (the 
intensities),  10  hidden  neurons  with  3  outputs  (r,  real  and  imaginary  part  of  m2). 

As  a  measure  of  goodness,  we  have  calculated  the  correlation  coefficient  defined  as  follows: 

where  xt  and  x2  are  two  observed  quantities  and  xj  is  the  mean  of  the  observed  value  for  the  i-type 
sample. 

In  Figs  3-6  we  collected  the  relevant  results  obtained  by  means  the  software  package  MATLAB 
implemented  on  a  Pentium  III  with  256  MB  of  RAM.  In  fig.4,  5  and  6  we  can  observe  that  after  1000 
training  cycles  (epochs)  the  retrieved  coefficients  (red  curves)  are  in  a  quite  good  agreement  with  the 
target  parameters  (blue  curves).  Even  if  the  better  agreement  regards  the  radius  that  gives  a 
correlation  coefficient  of  0.99574  (see  table  II),  we  can  assess  that  the  performance  of  the  network  is 
still  good  also  for  the  dielectric  constants.  It  is  useful  to  highlight  the  goodness  of  the  index  of 
performance  of  the  network  given  by  the  small  value  of  mean  squared  error  (see  Table  II). 
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Fig.  1  Pattern  of  a  set  of  three  non 
normalized  input  data 


Fig.  la  Pattern  of  the  mixed  and  normalized 
input  data  represented  in  fig  1. 
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Fig  2  -  Scheme  of  the  Neural  Network  implemented  with 
the  software  package  MATLAB 


N68171-01-M-5907 


Training -Blue  Goal-Black 


10 


Performance  is  0.00179615,  Goa!  is  1e-005 


14. 


Uscl  test 


Fig.  3  Development  of  the  rms  error  for  training 
and  test  data  set  during  neural-network  training 


Usc2  test 


Fig.  5  The  same  of  fig  3  but  for  the  real 
part  of  dilectric  constant. 


Fig.  4  Comparison  between  the  target  values  of 
the  radius  (blue  curve)  and  simulated  ones,  output 
of  the  neural  network  (red  curve). 
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Fig.  6  The  same  of  fig  3  but  for  the 
imaginary  part  of  dielectric  constant. 
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Efficient  Light  Scattering  Calculation  for  Aggregates  of  Large  Spheres 

It  is  well  known  that  a  possible  means  for  determining  shape  and  structure  information  about 
airborne  microparticles  is  to  investigate  the  angular  scattering  pattern  that  is  done  usually  in  only  one 
angular  dimension  (9). 

However,  as  the  elastic  scattered  angular  intensity  distribution  from  nonspherical  particles  has  an 
azimuthal  (9)  as  well  as  a  polar  (8)  dependence,  observation  of  the  two-dimensional  angular  optical 
scattering  (TAOS)  is  a  more  useful  tool  for  determining  the  shape  and  the  structure  of  microparticle 
aggregates.  Recently,  Holler  et  al.  (1)  have  devised  an  experimental  setup  that  is  able  to  illuminate  a 
single  scatterer  and  to  measure  its  TAOS  for  scattering  angles  that  go  from  forward  to  backward  direction. 
In  view  of  the  capability  of  this  new  experimental  device,  we  discuss  how  and  to  what  extent  the 
computational  effort  can  be  reduced  still  preserving  the  quantitative  features  of  the  signature  of  a  sample 
aggregate.  As  we  are  not  interested  to  use  orientational  averaging  procedure,  we  perform  all  the 
calculation  for  a  microparticle  aggregate  in  the  framework  of  the  E-scheme  that  prevents  us  to  obtain  the 
transition  matrix  that  requires  a  high  computational  effort. 

The  description  of  E-scheme,  the  results  of  the  calculation  and  the  comments,  also  presented  at 
the  Annual  Scientific  Conf.  On  Obscuration  ad  Aerosol  Res.  (Abardeen  Md.USA,  June  2001),  are  widely 
discussed  in  the  following  paper  that  will  be  published  in  Applied  Optics. 
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The  calculation  of  the  scattering  pattern  from  aggregates  of  spheres  through 
the  T-matrix  approach  yields  high  precision  results  but  at  an  high  computa¬ 
tional  cost  specially  when  the  aggregate  concerned  is  big  or  is  composed  of 
big-sized  spheres.  With  reference  to  a  specific  but  representative  aggregate, 
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still  preserving  the  qualitative  features  of  the  signature  of  the  aggregate 
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1.  Introduction 

The  scattering  properties  of  spherical  particles,  both  homogeneous  and  radially  non- 
homogeneous,  can  be  easily  and  accurately  calculated  through  Mie  theory1  or  one  of 
its  extensions2  over  a  wide  range  of  radii  and  wavelengths.  The  theoretical  predictions 
are  in  both  qualitative  and  quantitative  agreement  with  the  observations  provided  the 
experimental  setup  and  measured  parameters  faithfully  reproduce  the  assumed  theo¬ 
retical  model  (actual  sphericity  of  the  scatterer,  dielectric  constant,  scattering  angle, 
etc...).3 

The  situation  is  quite  different  when  the  scatterer  is  not  spherically  symmetric; 
in  that  case  the  calculation  of  the  scattering  properties  (scattering  amplitude,  scat¬ 
tering  and  extinction  cross  sections,  scattered  intensity)  with  an  accuracy  of  3  or  4 
significant  digits,  may  become  quite  a  cumbersome  task,  particularly  as  the  size  of 
the  scatterer  gets  bigger  and  bigger.  Performing  calculations  to  so  high  an  accuracy, 
however,  makes  sense  when  one  is  interested  in  the  properties  of  just  one  well  char¬ 
acterized  particle  but,  when  the  theoretical  predictions  are  to  be  compared  with  the 
experimental  data  from  an  assembly  of  scatterers,  such  a  computational  effort  may 
be  pointless.  Consider,  for  instance,  a  sample  made  of  a  dispersion  of  non  spherical 
particles:  as  a  rule,  the  component  particles  are  not  identical  nor  equioriented  so  that 


1 


r 


many  parameters  are  required  to  fully  characterize  the  dispersion.  In  view  of  the  un¬ 
avoidable  imprecision  in  the  determination  of  each  of  the  parameters  involved  in  the 
description  of  the  sample,  a  complete  characterization  of  the  scatterers  is  not  needed: 
a  model  reproducing  the  overall  shape  and  symmetry  of  the  scatterer  may  be  a  good 
one  to  perform  the  calculations.  Moreover,  the  effort  devoted  to  the  calculation  of 
the  scattering  properties  of  individual  particles  can  be  somewhat  reduced.  Of  course 
the  question  arises:  what  degree  of  precision  is  reasonable  to  pursue  in  the  calcu¬ 
lations  to  get  a  satisfactory  description  of  the  scattering  properties  of  the  particles 
and  what  is  the  most  suitable  computational  technique  to  achieve  this  precision?  An 
accuracy  within  e.g.  10%  could  be  enough  to  give  a  good  qualitative  description  of 
the  individual  scatterers,  to  be  integrated  in  the  overall  description  of  the  dispersion. 

Assuming  that  a  model  for  the  particles  has  been  chosen,  one  has  still  to  choose 
the  most  suitable  computational  technique  to  get  a  satisfactory  description  of  their 
scattering  properties,  as  several  computational  methods  imply  an  intrinsic  accuracy 
that  stems  from  the  assumptions  on  which  they  are  based.  In  this  respect,  methods 
based  on  series  expansions  are  able  to  give  results  whose  accuracy  can  be  graduated 
according  to  the  need:  in  fact,  we  recall  that  Mie  theory  itself  is  .just  a  series  ex¬ 
pansion  method.  Let  us  consider  also  the  problem  of  scattering  from  a  particle  that 
actually  is  or  can  be  modeled  as  an  aggregate  of  spherical  scatterers.  By  expanding 
the  field  in  a  series  of  vector  multipole  fields  one  is  able  to  get  the  transition  matrix 
and  the  scattering  amplitude  of  the  whole  aggregate  within  an  approximation  that 
depends  only  on  the  maximum  order  of  the  multipole  fields  that  are  included  into  the 
expansion.  This  is  a  useful  circumstance  because,  when  the  number  of  the  spheres  in 
the  aggregate  is  large  or  when  the  individual  size  of  the  component  spheres  is  big  the 
computational  effort  grows  more  and  more  heavy,  so  that,  determining  the  minimal 
acceptable  accuracy  of  the  results  becomes  an  important  issue. 

In  order  to  give  a  quantitative  content  to  the  above  considerations  we  devote 
this  paper  to  the  determination  of  the  minimal  acceptable  accuracy  of  scattering 
calculation  for  a  big  cluster  of  spheres  through  the  method  of  expansion  that  we 
mentioned  above.  To  this  end,  we  chose  to  deal  with  the  same  model  scatterer  that  has 
already  been  considered  by  Holler  et  al.4  We  stress,  however,  that  our  purpose  is  not  to 
compare  our  calculations  with  the  predictions  in  Ref.4  or  with  the  experimental  data; 
our  purpose  is  instead  to  investigate  through  the  specific  case  we  deal  with  how  and  to 
what  extent  computational  efforts  can  be  reduced  still  getting  a  reliable  qualitative 
description  of  the  scattering  pattern.  To  perform  our  task  a  two-step  calculation 
proved  to  be  more  viable  than  a  single-step  one.  Accordingly,  first  we  reduce  the 
geometrical  dimensions  of  the  scatterer  leaving  unchanged  everything  else;  doing  this, 
on  one  hand  drastically  reduces  the  difficulties  connected  with  the  computation  and 
the  time  required  to  perform  it,  and  on  the  other  hand  allows  us  to  work  out  the 
conditions  that  ensure  the  stability  of  the  results  from  a  qualitative  point  of  view.  In 
other  words  we  want  to  determine  the  conditions  under  which  the  calculated  scattering 
pattern  qualitatively  reproduce  the  exact  one  in  the  sense  that  it  reproduces  the 
scattering  signature  of  the  individual  particle.  In  the  second  step,  we  consider  the 
model  scatterer  with  its  full  dimensions  and  we  apply  the  already  tested  stability 
criteria  so  to  get  a  scattering  pattern  whose  validity,  at  least  from  a  qualitative  point 
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of  view,  is  reasonably  ensured. 

Section  2  gives  a  brief  description  of  the  theoretical  approach,  devoting  particular 
attention  to  the  aspects  regarding  the  convergence  and  the  stability  of  the  results. 
Section  3  is  devoted  to  the  presentation  and  discussion  of  the  results  and  finally  in 
section  4  we  draw  a  few  conclusive  remarks. 

2.  Scattering  by  clusters  of  spheres.  Theory 

In  order  to  fix  the  notation,  let  us  state  that  the  polarization  of  the  incident  and  of 
the  scattered  field  is  described  with  reference  to  two  pairs  of  mutually  orthogonal  unit 
vectors,  u/,,  and  Us,,,  with  rj  =  1,2,  such  that 

U/i  X  U/2  =  k/,  usi  x  US2  =  k s- 

In  this  paper  we  choose 

U/i  =  #/,  U/2  =  <fij 

and 

USl  =  #S,  «S2  =  <fisi 

where  #/,  tfij,  $s,  <Ps  are  unit  vectors  along  the  meridians  and  the  parallels  of  the 
unit  spherical  surface  with  center  within  the  scattering  particle  at  the  points  r  =  k/ 
and  r  =  ks,  respectively.5  Hereafter  we  denote  an  incident  field  polarized  along  U/r) 
as  E iv  and  the  corresponding  scattered  field  as  Esrj- 

Let  us  consider  a  cluster  of  spheres,  numbered  by  the  index  a,  of  radius  pa  and 
refractive  index  na ,  whose  centers  lie  at  Ra,  embedded  into  a  non-absorptive  medium 
of  refractive  index  n.  Neither  the  radii  nor  the  refractive  indexes  of  the  component 
spheres  need  to  be  equal  to  each  other.  The  electromagnetic  field  is  expanded  in  terms 
either  of  the  vector  multipole  fields  j{^,  with 

j£LV.  K)  =  f),  j£(r,  K)=tVx  jW(r,  K) 

that  are  regular  at  the  origin,  or  of  the  multipole  fields  H{^  that  are  identical  to  the  J 
multipole  fields  except  for  the  substitution  of  the  Hankel  functions  of  the  first  kind  hi 
in  place  of  the  Bessel  functions  ji  and  thus  satisfy  the  radiation  condition  at  infinity. 
In  the  preceding  equations  the  vector  functions  X/m  are  vector  spherical  harmonics8 
and  the  superscript  p  is  a  parity  index  that  distinguishes  the  magnetic  multipoles 
(p  =  1)  from  the  electric  ones  (p  =  2).  We  assume  that  all  the  fields  depend  on  time 
through  the  factor  exp[— iut\,  that  is  omitted  throughout,  and  define  the  propagation 
constant  in  vacuo  k  =  ui/c. 

The  incident  field  is  assumed  to  be  a  polarized  plane  wave  and  is  expanded  as 
E„(r)  =  £o£j£(r,n«0<L 

plm 
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where  Wj'jm  =  W$>( u,„k,)  with 

Wff(e,  k)  =  4^'-‘e  ■  Z“*(fi), 
where  we  define  the  transverse  harmonics 

z2(k)  =  Xlm((c),  zg(t)  =  Xj»(t)  x  fc. 

The  field  scattered  by  the  whole  cluster  is  written  as  the  superposition  of  the  fields 
scattered  by  the  single  spheres 

=  (1) 

a  plm 

where  ra  =  r  —  Ra,  whereas  the  field  within  the  a-th  sphere  is  expanded  as 

Er„0  =  S,^jS(ra,n.t)C£>m-  (2) 


The  amplitudes  A^lm  and  C^lm  in  Eqs.(l)  and  (2)  are  determined  by  the  custom¬ 
ary  boundary  conditions  across  the  surface  of  each  sphere  in  the  cluster,  and  their 
knowledge  totally  solves  the  problem  of  the  dependent  scattering  from  the  aggregate; 
in  fact,  for  each  value  of  a,  the  A^lm  describe  the  scattering  by  the  a-th  sphere, 
due  to  the  incident  plane  wave  and  to  the  waves  previously  scattered  by  all  the  other 
spheres  in  the  aggregate.  Anyway,  using  the  procedure  that  is  fully  described  e.g. 
in  Ref.  9,  the  amplitudes  A^lm  turn  out  to  be  the  solution  to  the  system  of  linear 
non-homogeneous  equations 

£  £  =  -<J L  (3) 

a'  p'/'m' 

where  we  define  the  shifted  amplitudes  of  the  incident  field 

\A)(p)  -  V  7(vp,)  W{p,)  ( 4 ) 

yVT)atlm  ~  JotlmOVm'VV VV 

p'l'm' 

The  quantities  Vm’  are  elements  of  the  matrix  that,  according  to  the  addition 
theorem  for  multipole  fields,6  translates  these  fields  from  the  origin  of  the  coordinates 
0  at  Ro  =  0  to  Rq,  i.e.  to  the  center  of  the  a-th  sphere,  and 

M^ln  w  =  (5) 

In  Eq.  (5)  the  quantities  R^}  and  are  the  elements  of  the  transition  matrix  for  the 
a-th  sphere  and,  except  for  a  sign,  coincide  with  the  customary  Mie  coefficients  6/  and 
ai ,  respectively;  the  quantities  whose  explicit  expression  is  given  in  Ref.7, 

come  from  the  addition  theorem  of  Ref.  6  take  account  of  the  multiple  scattering 
processes  that  occur  among  the  spheres  in  the  aggregate. 
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A .  T-scheme 

In  principle,  the  amplitudes  contain  all  the  information  on  the  scattering 

properties  of  the  cluster  concerned.  In  fact,  we  will  now  discuss  two  schemes  for  the 
extraction  of  the  relevant  information  from  the  solution  to  system  (3),  the  choice 
being  dictated  by  the  kind  of  information  one  is  interested  in.  The  first  scheme,  that 
we  summarize  in  the  present  subsection,  is  the  most  complete  as  it  produces  the  T- 
matrix  of  the  scatterer.  To  this  end,  letting  M  be  the  matrix  whose  elements  are  given 
by  Eq.(5),  we  recall  that  the  formal  solution  to  system  (3), 

=  -£  [M-t <•> 

p'/'m' 

resembles  the  definition  of  the  T-matrix.10  Nevertheless,  M_1  does  not  coincide  with 
the  T-matrix  of  the  aggregate  because,  according  to  Ref.  10,  the  transition  matrix 
relates  the  multipole  amplitudes  of  the  incident  field  to  those  of  the  field  scattered  by 
the  whole  object,  whereas  Eq.  (6)  relates  the  amplitudes  of  the  incident  field  to  those 
of  the  fields  scattered  by  each  sphere  in  the  aggregate.  Nevertheless,  the  addition 
theorem  of  Ref.6  allows  us  to  write  the  scattered  field  in  terms  of  multipole  fields 
with  origin  at  O  as 

e*,  =  £>.££  £ 

p/m  a*  p'/'m' 

and  consequently  to  write  the  amplitudes  of  the  field  scattered  by  the  whole  aggregate 
as 

a(p)  -  v  V  7{ppt)  A{pt)  (7) 

a1  p'/'m' 


where  the  quantities  whose  explicit  expression  is  given  in  Ref.  7,  are  the 

elements  of  the  matrix  that  translates  the  origin  of  the  H-multipole  fields  from  Ra/ 
to  the  origin  of  the  coordinates  at  O. 

At  this  stage  the  expression  for  the  amplitudes  given  by  Eq.  (6)  can  be 

substituted  into  Eq.  (7)  that  consequently  reads 


=  £4 


o>f') 

/m/W  KK/y/'m'> 


p' I'm ' 


where  the  elements  of  the  true  T-matrix  of  the  aggregate  are 


o(pp')  _ 


=  -£££ 

aa>  qLMq'L'M ‘ 


<t(p?) 

^ QlmaLM  [ 


i-i]  (kj') 

J  aLMa'L'M' 


•Ja'L'M'OV  m' 


(8) 


The  elements  account  for  all  the  scattering  properties  of  the  cluster.  In  fact, 

in  the  far  zone  the  scattered  field  assumes  the  asymptotic  form 


Esrj  =  Eq 


exp  (inkr) 


f^(ks,k/) 


Eq  exp  (inkr) 
nk  r 


£h)'wz  !5(k)4SL 
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where  is  the  scattering  amplitude8  and  Eq.  (7)  has  been  used.  Therefore  the  com¬ 
ponent  of  the  scattered  field  that  is  polarized  along  u sn  reads 

expifcr 

EJSr}'  *  ttSVj  —  Jr)T}/ 

r 

with 


fi m'  — 


4irnk 


Y''  v™'  w(p)*  w(p') 

/  ,  /  ,  VV ST)lm°lml’mlVV  In'Vm' 


(9) 


plm  p'Vm’ 


=  w£\uSnM  As  a  consequence,  the  scattered  intensity,  both  co¬ 
polarized  and  cross-polarized,  turns  out  to  be 

Hereafter  we  will  refer  to  the  procedure  outlined  above  as  the  T-scheme. 


B.  E-scheme 

We  now  come  to  describe  an  alternative  scheme,  the  El-scheme,  that,  in  some  in¬ 
stances,  has  useful  computational  advantages.  We  start  remarking  that,  when  Eq.  (8) 
is  substituted  into  Eq.  (9)  to  get  the  components  of  the  scattering  amplitude,  the 
resulting  expression  contains  the  terms 


E 

p'/'m' 


'jWv')  li/b7') 
Ja'L'WQVm'™  In' l'ml 


that  can  be  somewhat  simplified  as  follows.  Since  r  =  Ra  +  ra, 

exp[zk/  •  r]  =  exp[ik/  •  RQ]  exp[zkj  •  ra]  . 

Therefore,  by  expanding  u/^exp^k/  •  r]  and  Uj^exp^'k/  •  ra\  we  get 

wSL = »p[*k/  • 

Comparison  of  the  latter  equation  with  Eq.  (4)  yields 

E  = «piik'  ■  rj<L-  <10) 

p* Vm 1 

Substitution  into  Eq.(10)  of  Us,,  in  place  of  U/,,  and  kg  in  place  of  k/  and  complex 
conjugation  then  yields 

E  =  exp[-iks  • 

•pWm* 

Now,  on  account  that  the  refractive  index  n  is  real,  we  get11 

_  7 (p'p) 

'JalmOl'm'  ~  <^0 I'm' aim' 
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so  that 


E  wsrX'^£.lm  =  exP[~*ks  ■  RJ<t.  (11) 

p,lfml 

Using  Eqs.(10)  and  (11),  fm*  takes  on  the  final  form 

fnn'  =  ^  E  E  E  Ws%  exP[— *ks  •  R*] 

K  pLMp'UM '  aa' 

X  l]^LMa'L'M'  exP[*k/  •  Rq'] 

We  stress  that  unlike  the  approximate  scheme  of  Xu12  to  get  the  expression  of  the 
scatterd  field  in  the  far  zone,  the  El-scheme  is  exact  as  it  implies  no  further  approxi¬ 
mation  beyond  the  customary  truncation  of  the  multipole  expansions.  We  tested  the 
reliability  of  this  scheme  by  using  it  extensively  to  compare  our  theoretical  predic¬ 
tions13  with  the  experimental  data  of  Schuerman  and  Wang  from  well  characterized 
clusters  in  fixed  orientation.14,15 

C.  Convergence 

According  to  T-scheme,  calculation  of  the  transition  matrix  of  a  cluster  requires 
inverting  the  matrix  M  whose  order  is,  in  principle,  infinite.  Naturally,  system  (3)  must 
be  truncated  to  some  finite  order  by  including  in  Eq.(8)  terms  up  to  order  Lm,  that  is 
the  maximum  value  both  for  L  and  V .  chosen  so  as  to  ensure  the  required  accuracy  of 
the  transition  matrix  elements.  For  a  cluster  of  N  spheres  this  implies  to  solve  a  system 
of  order  Dm  =  2NLm(Lm  +  2)  that  may  grow  rather  big.  Actually,  the  inversion  of 
the  matrix  M  is  responsible  for  most  of  the  time  required  for  the  calculation;  this 
time  scales,  in  fact,  as  Dm-  According  to  the  considerations  in  Section  1,  Lm  must 
be  so  chosen  that  using  an  increased  Lm  the  pattern  of  the  scattered  intensity  no 
more  undergoes  any  qualitative  changes.  Choosing  the  least  adequate  Lm  is  thus 
decisive:  on  account  of  the  definition  of  Dm,  computation  time  increases  as  L6M  and 
storage  as  L4M .  Therefore,  using  too  big  an  Lm  would  require  a  computer  time  that  is 
definitely  longer  than  a  real-time  calculation  that  may  be  required  e.g.  for  medical  or 
meteorological  applications.  It  may  be  worth  remarking  that  Lm  could  take  a  different 
value  for  each  of  the  spherical  components  of  the  cluster  we  deal  with:  in  other  words, 
there  is  no  reason  why  l  in  Eqs.(l)  and  (2)  should  not  assume  a  maximum  value  that 
depends  on  a.  As  a  noticeable  consequence,  Dm  could  be  somewhat  reduced.  This 
is  in  connection  with  the  obvious  remark  that  the  spherical  components,  according 
to  their  own  scattering  power,  may  give  considerably  different  contributions  to  the 
scattering  from  the  cluster  as  a  whole.  Such  a  refinement  in  assigning  Lm,  however, 
is  of  no  use  when  the  scattering  powers  of  the  individual  components  are  similar  and 
when  one  deals  with  compact  structures.  Since  in  Sect.3  we  will  deal  just  with  a 
cluster  of  this  kind,  we  will  assign  only  one  value  of  Lm  to  be  applied  to  each  of  the 
spheres  composing  the  cluster. 

At  this  stage,  it  may  be  useful  to  make  a  few  quantitative  considerations  about 
choosing  Lm  for  exact  convergence,  i.e.  convergence  to  fa  3-4  significant  digits.  For 
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a  dielectric  sphere  with  radius  ba  we  have  to  choose  Lm  >  kba  +  ma;  when  we  deal 
with  a  cluster  of  spheres,  we  get  the  exact  description  of  the  scattering  problem 
considering  Lm  >  kbc  +  mc,  where  bc  is  the  radius  of  the  smallest  sphere  including 
the  whole  aggregate.  The  numbers  ma  and  mc  are  positive  integers  that  depend  on 
the  refractive  indexes;  in  general,  if  kba  and  kbc  are  large  ma  <C  kba  and  mc  <C  kbc. 
Nevertheless,  in  the  next  section,  we  will  show  that,  in  order  to  give  a  good  qualitative 
description  of  the  scattering  properties  of  the  aggregate  we  consider  it  is  sufficient  to 
use  a  value  for  Lm  that  would  rather  be  appropriate  to  the  exact  description  of  the 
scattering  from  a  single  sphere  of  the  aggregate.  The  validity  of  our  choice  for  Lm  will 
be  checked  by  looking  wether  or  not  the  scattering  properties  of  the  whole  aggregate, 
thus  calculated,  will  still  keep  their  qualitative  features  with  increasing  value  of  Lm- 

Finally,  let  us  stress  that  the  T-scheme  and  the  El-scheme  are  not  the  only  pos¬ 
sible  ones  to  describe  scattering  from  aggregated  spheres.  In  fact,  schemes  have  been 
devised  that  do  not  require  direct  inversion  of  matrix  M  as  well  as  schemes  in  which 
the  averaging  and  inversion  are  integrated  in  a  single  computational  procedure.16,17 
As  all  these  alternative  schemes  are  aimed  at  minimizing  the  computational  effort, 
they  should  be  carefully  considered  when  the  main  interest  is  the  exact  solution  of 
the  scattering  problem. 

Using  the  El-scheme  slightly  reduces  the  amount  of  calculation  for  a  given  Lm 
with  respect  to  the  T-scheme.  In  fact,  the  latter  requires  the  additional  calculation 
of  the  sums  in  the  left  hand  side  of  Eqs.(10)  and  (11)  that,  in  turn,  is  more  or  less 
onerous  according  to  the  value  one  assumes  for  Im ,  the  maximum  l1  to  be  included 
to  get  an  accurate  representation  of  the  scattered  field.  El-scheme,  however,  prevents 
us  to  obtain  the  transition  matrix  as  defined  by  Eq.  (8)  and  therefore  to  use  the 
orientational  averaging  procedure,  widely  described  elsewhere:18  the  true  transition 
matrix,  indeed,  cannot  be  reminiscent  of  information  either  on  the  excitation  (through 
the  factors  exp[ikj  •  RQ])  or  on  observation  (through  the  factors  exp[— iks  •  Ra<]).  It 
is  also  to  be  borne  in  mind  that  once  the  T-scheme  has  been  applied,  another  chunk 
of  calculation  is  needed  to  perform  the  orientational  averages.  In  conclusion  the  El- 
scheme  is  appropriate  when  orientational  averages  are  not  required. 

3.  Results  and  discussion 

Although  the  cluster  that  we  study  in  this  paper  is  the  same  13-spheres  clusters 
considered  by  Holler  et  al.,  we  stress  again  that  our  purpose  is  quite  different,  as  we 
explained  in  Sect.l.  The  cluster  of  interest  is  composed  of  polystirene  latex  spheres 
with  size  parameter  xa  =  nkb,  =  13.36,  arranged  in  a  closely  packed  structure  with 
a  roundish  geometry  so  that  the  overall  size  parameter  of  the  whole  aggregate  turns 
out  to  be  xc  =  nkbc  =  40.  The  coordinates  of  the  centers  of  the  spheres  are  listed  in 
Table  1;  n  =  1.  According  to  the  considerations  at  the  end  of  Sect.  2  convergence  to 
3  decimal  places  would  require  Lm  >40.  In  Ref.  4,  the  calculations  were  performed 
through  the  recursive  algorithm  by  Chew19  that  shares  in  common  with  our  approach 
the  fact  that  the  only  approximation  lies  in  the  choice  of  the  value  of  Lm  at  which  the 
multipole  expansions  of  the  electromagnetic  fields  are  stopped.  Nevertheless  it  appears 
that  a  generalized  use  of  such  a  perturbative-iterative  approach  would  require  some 
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caution  and  artwork  in  its  implementation.  For  instance,  the  actual  achievement  of 
convergence  for  a  given  choice  of  Lm  should  be  carefully  verified.  On  the  contrary, 
the  direct  approach  that  we  use,  i.e.  the  El-scheme,  appears  to  be  straightforward  and 
leaves  to  the  researcher  full  control  of  the  accuracy  that  can  be  achieved.  ^ 

In  the  following  discussion  the  wavevector  of  the  incident  plane  wave  k/  has 
(pj  =  0°  and  i)i  =  90°,  or  i?/  =  45°.  The  scattered  intensity  is  calculated  for  -90°  < 
•ds  <  90°  and  0°  <  <ps  <  180°  and  is  represented  through  contour  plots.  In  a  few  cases 
are  also  reported  sections  of  these  contour  plots  for  ds  =  30°  and  ds  —  90  .  Choosing 
di  =  90°  corresponds  to  incident  field  viewing  the  scatterer  as  a  highly  symmetric 
object. 

With  the  aim  of  finding  a  qualitative  stability  criterion  we  consider  a  cluster  of 
spheres  with  the  same  geometry  as  the  one  we  chose  but  scale  down  the  dimensions 
of  the  component  spheres  to  the  20  %.  Since  we  assume  that  the  incident  field  has  a 
wavelength  of  0.532  /um,  the  size  parameter  will  be  reduced  to  xa  =  2.672  for  the  single 
component  sphere  and  to  xc  =  8  for  the  whole  cluster.  We  make  an  angular  analysis 
of  the  scattered  intensity  in  the  forward  scattering  zone.  In  Fig.  1  (a),  (b)  and  (c)  we 
report  the  contour  plots,  calculated  using  the  E-scheme,  of  the  differential  scattering 
cross  section  =  |/,w|2  for  lM  =  3,  5,  8,  respectively.  The  corresponding  plots  for 
Xw  show  the  same  trend  that  we  are  going  to  discuss  and  are  therefore  not  reported. 
Comparison  of  the  plots  (a),  (b)  and  (c)  shows  that  the  stability  is  already  reached 
for  Lm  =  5;  in  fact,  the  plots  for  LM  =  5  [Fig.l  (b)]  and  LM  =  8  [Fig.l  (c)]  do  not 
show  any  important  difference.  Note  that  LM  =  5  is,  in  principle,  a  little  larger  than 
barely  sufficient  to  get  the  convergency  for  the  component  spheres.  This  result  has 
been  further  investigated  by  reporting  in  Fig.  2  (a)  the  differential  scattering  cross 
section  for  dj  =  90°  and  <pi  =  0°  for  Lm  =  2,  3,  4,  5;  in  Fig. 2  (a)  ds  =  90  and 
-90°  <  (ps  <  90°;  in  Fig.  2(b)  ds  =  30°  and  -90°  <  <  90°.  The  plots  in  both 

Fig.2  (a)  and  Fig.2  (b)  show  that  the  changes  undergone  by  the  intensity  profiles  for 
Lm  >  4  are  quite  negligible.  We  do  not  report  the  profiles  for  Lm  >  5  because,  on 
the  scale  of  the  figure  they  are  perfectly  superposed  to  the  profile  for  Lm  =  5-  We  are 
thus  led  to  conclude  that  the  contour  plot  in  Fig.l  (b),  i.e.  for  Lm  =  5,  gives  a  good 
description  of  the  actual  scattered  intensity  and  that  Lm  =  4  is  sufficient  to  get  an 
acceptable  qualitative  description. 

As  regards  the  behavior  of  the  T-scheme,  we  found,  as  expected,  that  Im  >  Lm 
is  required  to  get  full  convergence  in  the  sums  in  Eq.(10)  and  (11)  so  as  to  reproduce 
the  results  obtained  in  the  E-scheme  for  the  same  value  of  Lm-  In  practice,  we  report 
in  Fig. 3  the  intensity  profiles  of  Fig.2  (a),  calculated  with  Lm  —  4  both  in  the  E  and 
T-scheme  and,  in  the  latter  case,  for  increasing  values  of  lM-  The  curves  in  Fig.3  show 
at  once  that,  although  Im  =  8  is  quite  acceptable,  and  that  Im  =  10  is  even  better, 
the  results  of  the  E  and  T-scheme  become  indistinguishable  on  the  scale  of  the  figure 
for  Im  =  12.  Such  a  perfect  convergence  is  unnecessary,  however,  as  it  would  yield 
the  perfect  reproduction  of  the  results  of  the  E-scheme  for  Lm  =  4,  that  are  in  their 
nature  qualitative  results.  Thus,  one  could  choose  for  Im  the  least  value  that  would  not 
disrupt  the  qualitative  behavior  of  the  scattering  pattern:  according  to  Fig.3,  Im  =  8 
is  a  fair  choice.  This  is  yet  a  large  value  so  that  there  is  a  cost,  in  particular  when 
averaging  over  the  orientation  of  the  scatterer  is  in  order.  Anyway,  we  can  conclude 
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that  using  the  E-scheme,  a  good  stability  criterion  requires  using  a  value  of  Lm  just 
a  little  bigger  than  necessary  to  get  the  convergence  for  the  component  spheres. 

In  order  to  test  the  stability  criterion  just  found,  we  calculate,  using  the  E-scheme, 
the  scattering  pattern  of  the  cluster  considered  with  its  real  dimensions:  in  this  respect 
we  recall  that  now  xs  =  13.36.  The  results  are  reported  for  dj  —  90°  and  <pj  =  0° 
in  Fig.4  for  Lm  =  13  in  (a)  and  for  Lm  =  15  in  (b);  the  pattern  for  dj  =  45°  and 
<pj  =  0°  is  reported  in  Fig. 5  for  Lm  =  13  in  (a)  and  for  Lm  =  15  in  (b).  Even  a 
cursory  examination  of  the  patterns  in  Figs. 4  and  5  shows  that  the  stability  criterion 
is  effective  in  yielding  patterns  that  give  an  acceptable  qualitative  description  of  the 
scattering  properties  of  the  cluster  at  hand.  In  fact,  when  Lm  is  increased  from  13  to 
15  the  changes  undergone  by  the  patterns  are  small  on  the  scale  of  the  figure.  Even  in 
this  case,  however,  the  numerical  output  shows  that  an  accuracy  to  3  decimal  places 
has  not  yet  been  achieved. 

4.  Conclusions 

At  first  glance,  the  stability  criterion  that  we  found  in  Sect.  3  may  be  rather  surprising. 
In  fact,  the  calculations  in  Sect.  3  support  the  conclusion  that  by  chosing  an  Im  barely 
sufficient  to  ensure  the  convergence  of  the  field  scattered  by  the  component  spheres, 
one  also  achieves  an  acceptable  stability  of  the  scattering  pattern  from  the  whole 
cluster.  Nevertheless,  let  us  recall  that  the  pattern  from  the  cluster  differ  from  that  of 
the  independent  spheres  because  of  the  coupling  effected  by  the  multipole  translation 
matrix  elements  Eq.  (5).  These  quantities,  as  appears  from  their  explicit 

expression  in  Ref.  7,  heavily  depend  on  the  distance  between  the  centers  of  the  spheres 
and  are  bound  to  be  small,  even  for  contacting  spheres,  when  the  radius  of  the  latter 
is  large.  Of  course,  the  translation  matrix  elements  are  even  smaller  for  noncontacting 
spheres.  As  a  result,  in  spite  of  the  coupling,  the  scattering  pattern  from  a  cluster  of 
big  spheres  is  governed  by  the  scattering  power  of  the  spherical  components,  that  is 
given  by  the  Mie  coefficients  R^.  Of  course,  the  criterion  that  we  found  in  Sect.  3,  may 
not  be  strictly  valid  for  a  cluster  whose  spherical  components  have  strongly  different 
scatting  power.  In  such  a  case  a  more  detailed  analysis  is  in  order  by  chosing,  as  we 
suggested  in  Sect.  2,  a  different  Im  for  each  sphere. 

Anyway,  the  calculations  that  we  presented  in  Sect.  3  support  our  considerations 
in  Sect.  1  on  the  accuracy  that  is  reasonable  to  pursue  when  calculating  the  optical 
properties  of  model  particles.  In  fact,  the  particles  that  are  produced  in  actual  ex¬ 
periments  are  often  not  well  characterized  whereas  the  mathematical  techniques  may 
give  results  that  are  strongly  dependent  on  the  exact  geometry,  size  and  orientation  of 
the  particles  concerned.  This  is  just  the  problem  that  has  been  faced  by  Mishchenko 
and  Mackowski  in  comparing  their  theoretical  predictions  with  the  experimental  data 
from  bispheres.20’21  These  authors  state  that  they  were  forced  to  repeat  their  calcula¬ 
tions  several  times  changing  the  size  of  the  particles  in  order  to  get  results  reasonably 
resembling  the  experimental  findings.  It  is  evident  that  the  computational  effort  in 
the  case  of  bispheres  is  much  lighter  than  in  the  case  of  the  13-sphere  cluster  we 
deal  with  in  this  paper.22  Therefore,  using  the  stability  criterion  we  established  in 
Sect  3,  possibly  in  the  framework  of  the  E-scheme,  will  save  a  lot  of  calculations, 
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still  reproducing  the  main  features  of  the  signature  of  the  particle  concerned.  This, 
in  turn,  would  make  viable  repeating  of  the  calculations  in  order  to  get  the  patterns 
that  better  resemble  the  experimental  data. 
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List  of  Figure  Captions 

Fig.  1.  Contour  plot  of  the  differential  scattering  cross  section  (log  scale)  of  the  ag¬ 
gregate  reduced  to  20  %  of  its  real  size,  calculated  in  the  E-scheme.  The  incident  field 
has  $]  =  90°  and  <pi  =  0°.  In  (a)  Lm  =  3,  in  (b)  Lm  =  5  and  in  (c)  Lm  =  8. 

Fig.  2.  Differential  scattering  cross  section  for  the  cluster  reduced  to  20  %  of  its  real 
size,  calculated  in  the  E-scheme  for  Lm  =  2  (dotted),  3  (dot-dashed),  4  (bulleted) 
and  5  (solid),  for  i?/  =  90°  and  pi  =  0°.  In  (a)  i?s  =  90°;  in  (b)  i?s  =  30°. 

Fig.  3.  Comparison  of  the  differential  scattering  cross  section  for  the  cluster  reduced 
to  20  %  of  its  real  size  calculated  in  the  E-scheme  for  Lm  =  4  (heavy  solid  line)  with 
that  yielded  by  the  T-scheme  with  Lm  —  4  and  Im  =  4  (light  solid),  6  (dot-dashed), 
8  (dotted),  10  (bulleted),  12  (dashed),  di  =  90°,  <pi  =  0°,  and  =  90°.  Note  that 
the  heavy  solid  curve  and  the  dashed  curve  are  almost  perfectly  superposed  on  the 
scale  of  the  figure. 

Fig.  4.  Contour  plot  of  the  differential  scattering  cross  section  for  the  cluster  con¬ 
sidered  with  its  real  size  calculated  in  the  E-scheme  for  i?/  =  90°  and  <pi  =  0°  with 
Lm  =  13  in  (a)  and  Lm  =  15  in  (b). 

Fig.  5.  Contour  plot  of  the  differential  scattering  cross  section  for  the  cluster  con¬ 
sidered  with  its  real  size  calculated  in  the  E-scheme  for  i?/  =  45°  and  <pi  =  0°  with 
Lm  =  13  in  (a)  and  Lm  =  15  in  (b). 
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Table  1.  Coordinates  of  the  centers  of  the  component  spheres  in  fim.  The 
radius  and  refractive  index  of  each  sphere  are  1.145 /rm  and  1.59,  respectively. 


X 

y 

z 

1 

0.00000 

0.0000000 

0.0000000 

2 

0.00000 

1.9250000 

-1.1890000 

3 

0.00000 

-1.9250000 

-1.1890000 

4 

0.00000 

-1.9250000 

1.1890000 

5 

0.00000 

1.9250000 

1.1890000 

6 

1.92500 

-1.1890000 

0.0000000 

7 

-1.92500 

1.1890000 

0.0000000 

8 

-1.92500 

-1.1890000 

0.0000000 

9 

1.92500 

1.1890000 

0.0000000 

10 

-1.18900 

0.0000000 

1.9250000 

11 

1.18900 

0.0000000 

1.9250000 

12 

1.18900 

0.0000000 

-1.9250000 

13 

-1.18900 

0.0000000 

-1.9250000 

15 


